library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(base)
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
##
## first, last
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
dataPath <- "/Users/bianca/Documents/U of Chicago/Classes/Time Series"
contracts_volume <- read.csv(paste(dataPath,'Contracts_Volume.csv',sep = '/'), header=TRUE)
contracts_class <- read.csv(paste(dataPath,'Contracts_Classification.csv',sep = '/'), header=TRUE)
contracts_volume$Total.Volume <- as.numeric(contracts_volume$Total.Volume)
contracts_volume$Electronic.Volume <- as.numeric(gsub(",", "", contracts_volume$Electronic.Volume))
contracts_volume$Floor.Volume <- contracts_volume$Total.Volume - contracts_volume$Electronic.Volume
joined_data <- contracts_volume %>%
left_join(contracts_class, by = c("Commodity.Indicator" = "Commodity.Code"))
cme_data <- joined_data
cme_data$Division <- "CME"
cme_data$Date <- as.Date(cme_data$Date, format= "%m/%d/%Y")
imm_data <- filter(joined_data, Division == "IMM")
imm_data$Date <- as.Date(imm_data$Date, format= "%m/%d/%Y")
iom_data <- filter(joined_data, Division == "IOM")
iom_data$Date <- as.Date(iom_data$Date, format= "%m/%d/%Y")
cme_agg <- cme_data %>%
group_by(Commodity.Indicator, Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
imm_agg <- imm_data %>%
group_by(Commodity.Indicator, Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
iom_agg <- iom_data %>%
group_by(Commodity.Indicator, Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
print(cme_agg)
## # A tibble: 21,048 × 5
## # Groups: Commodity.Indicator [559]
## Commodity.Indicator Date Total_Electronic.Volume Total_Total.Volume
## <chr> <date> <dbl> <dbl>
## 1 16E 2011-02-01 11 11
## 2 16E 2011-03-01 24 24
## 3 16E 2011-04-01 85 85
## 4 16E 2011-05-01 63 63
## 5 16E 2011-06-01 70 70
## 6 16E 2011-07-01 44 44
## 7 16E 2011-08-01 8 8
## 8 16E 2011-09-01 12 12
## 9 16E 2011-10-01 5 5
## 10 1A 2007-05-01 0 116
## # ℹ 21,038 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>
print(imm_agg)
## # A tibble: 5,723 × 5
## # Groups: Commodity.Indicator [78]
## Commodity.Indicator Date Total_Electronic.Volume Total_Total.Volume
## <chr> <date> <dbl> <dbl>
## 1 16E 2011-02-01 11 11
## 2 16E 2011-03-01 24 24
## 3 16E 2011-04-01 85 85
## 4 16E 2011-05-01 63 63
## 5 16E 2011-06-01 70 70
## 6 16E 2011-07-01 44 44
## 7 16E 2011-08-01 8 8
## 8 16E 2011-09-01 12 12
## 9 16E 2011-10-01 5 5
## 10 36E 2011-03-01 27 27
## # ℹ 5,713 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>
print(iom_agg)
## # A tibble: 16,399 × 5
## # Groups: Commodity.Indicator [475]
## Commodity.Indicator Date Total_Electronic.Volume Total_Total.Volume
## <chr> <date> <dbl> <dbl>
## 1 1A 2007-05-01 0 116
## 2 1A 2008-01-01 0 25
## 3 1A 2009-04-01 0 14
## 4 1A 2011-03-01 76 426
## 5 1A 2011-06-01 26 26
## 6 1A 2011-07-01 0 220
## 7 1A 2011-08-01 11 11
## 8 1A 2011-09-01 4 4
## 9 1A 2011-11-01 146 146
## 10 1A 2011-12-01 25 28
## # ℹ 16,389 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>
cme_agg2 <- cme_data %>%
group_by(Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
imm_agg2 <- imm_data %>%
group_by(Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
iom_agg2 <- iom_data %>%
group_by(Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
imputed_cme <- read.csv(paste(dataPath,'imputed_cme.csv',sep = '/'), header=TRUE)
names(imputed_cme) <- c("Date", "Price")
imputed_cme$Date <- as.Date(imputed_cme$Date)
imputed_imm <- read.csv(paste(dataPath,'imputed_imm.csv',sep = '/'), header=TRUE)
names(imputed_imm) <- c("Date", "Price")
imputed_imm$Date <- as.Date(imputed_imm$Date)
imputed_iom <- read.csv(paste(dataPath,'imputed_iom.csv',sep = '/'), header=TRUE)
names(imputed_iom) <- c("Date", "Price")
imputed_iom$Date <- as.Date(imputed_iom$Date)
cme_summary <- left_join(imputed_cme, cme_agg2, by = "Date")
imm_summary <- left_join(imputed_imm, imm_agg2, by = "Date")
iom_summary <- left_join(imputed_iom, iom_agg2, by = "Date")
cme_numeric_data <- cme_summary %>%
select_if(is.numeric)
cme_cor_matrix <- cor(cme_numeric_data, use = "complete.obs", method = "pearson")
cme_cor <- as.data.frame(cme_cor_matrix)
imm_numeric_data <- imm_summary %>%
select_if(is.numeric)
imm_cor_matrix <- cor(imm_numeric_data, use = "complete.obs", method = "pearson")
imm_cor <- as.data.frame(imm_cor_matrix)
iom_numeric_data <- iom_summary %>%
select_if(is.numeric)
iom_cor_matrix <- cor(iom_numeric_data, use = "complete.obs", method = "pearson")
iom_cor <- as.data.frame(iom_cor_matrix)
print(cme_cor)
## Price Total_Electronic.Volume Total_Total.Volume
## Price 1.00000000 0.6530748 0.6802748
## Total_Electronic.Volume 0.65307484 1.0000000 0.9736634
## Total_Total.Volume 0.68027480 0.9736634 1.0000000
## Total_Floor.Volume -0.06265014 -0.3726483 -0.1512653
## Total_Floor.Volume
## Price -0.06265014
## Total_Electronic.Volume -0.37264828
## Total_Total.Volume -0.15126533
## Total_Floor.Volume 1.00000000
print(imm_cor)
## Price Total_Electronic.Volume Total_Total.Volume
## Price 1.0000000 0.3287382 0.4764459
## Total_Electronic.Volume 0.3287382 1.0000000 0.9465020
## Total_Total.Volume 0.4764459 0.9465020 1.0000000
## Total_Floor.Volume 0.2972396 -0.4706299 -0.1607259
## Total_Floor.Volume
## Price 0.2972396
## Total_Electronic.Volume -0.4706299
## Total_Total.Volume -0.1607259
## Total_Floor.Volume 1.0000000
print(iom_cor)
## Price Total_Electronic.Volume Total_Total.Volume
## Price 1.00000000 -0.02515719 0.05596529
## Total_Electronic.Volume -0.02515719 1.00000000 0.98393410
## Total_Total.Volume 0.05596529 0.98393410 1.00000000
## Total_Floor.Volume 0.43943993 -0.29485138 -0.11951925
## Total_Floor.Volume
## Price 0.4394399
## Total_Electronic.Volume -0.2948514
## Total_Total.Volume -0.1195193
## Total_Floor.Volume 1.0000000
It looks like Floor Volume is not a good indicator for CME class seats prices. Both Electronic Volume and Total Volume are not good indicators for IOM class seats prices. We can also see that Electronic Volume and Total Volume are highly correlated in all seat classes, so only one of those will be used as predictor.
train_data_cme <- subset(cme_summary, Date < "2013-01-01")
test_data_cme <- subset(cme_summary, Date >= "2013-01-01")
train_data_imm <- subset(imm_summary, Date < "2013-01-01")
test_data_imm <- subset(imm_summary, Date >= "2013-01-01")
train_data_iom <- subset(iom_summary, Date < "2013-01-01")
test_data_iom <- subset(iom_summary, Date >= "2013-01-01")
price_cme <- ts(train_data_cme$Price, frequency = 12)
price_imm <- ts(train_data_imm$Price, frequency = 12)
price_iom <- ts(train_data_iom$Price, frequency = 12)
lm_cme <- lm(Price ~ Total_Total.Volume, data = train_data_cme)
lm_forecast_cme <- predict(lm_cme, newdata = test_data_cme, n.ahead = 12)
summary(lm_cme)
##
## Call:
## lm(formula = Price ~ Total_Total.Volume, data = train_data_cme)
##
## Residuals:
## Min 1Q Median 3Q Max
## -419623 -112256 -7097 82349 791934
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.166e+05 4.127e+04 2.826 0.0054 **
## Total_Total.Volume 2.945e-03 2.371e-04 12.420 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 197600 on 142 degrees of freedom
## Multiple R-squared: 0.5207, Adjusted R-squared: 0.5173
## F-statistic: 154.3 on 1 and 142 DF, p-value: < 2.2e-16
lm_imm <- lm(Price ~ Total_Total.Volume, data = train_data_imm)
lm_forecast_imm <- predict(lm_imm, newdata = test_data_imm, n.ahead = 12)
lm_iom <- lm(Price ~ Total_Floor.Volume, data = train_data_iom)
lm_forecast_iom <- predict(lm_iom, newdata = test_data_iom, n.ahead = 12)
For CME linear regression, Electronic Volume and Total Volume are correlated so only Total Volume is used. Floor Volume doesn’t have high correltion with price so it’s not used as predictor. For IMM linear regression, Electronic Volume and Floor Volume is picked instead so there’s no overlap between Total Volume and Floor Volume as predictors. For IOM, only Floor Volume is highly correlated with price so it’s the only predictor.
arma_model_cme <- auto.arima(price_cme, xreg = train_data_cme$Total_Total.Volume)
arma_forecast_cme <- forecast(arma_model_cme, xreg = test_data_cme$Total_Total.Volume, h=12)
arma_model_imm <- auto.arima(price_imm, xreg = train_data_imm$Total_Total.Volume)
arma_forecast_imm <- forecast(arma_model_imm, xreg = test_data_imm$Total_Total.Volume, h=12)
arma_model_iom <- auto.arima(price_iom, xreg = train_data_iom$Total_Total.Volume)
arma_forecast_iom <- forecast(arma_model_iom, xreg = test_data_iom$Total_Total.Volume, h=12)
library(forecast)
hw_model_cme <- hw(price_cme, xreg = train_data_cme$Total_Total.Volume)
hw_forecast_cme <- forecast(hw_model_cme, xreg = test_data_cme$Total_Total.Volume, h = 12)
hw_model_imm <- hw(price_imm,xreg = train_data_imm$Total_Total.Volume)
hw_forecast_imm <- forecast(hw_model_imm, xreg = test_data_imm$Total_Total.Volume, h = 12)
hw_model_iom <- hw(price_iom, xreg = train_data_iom$Total_Total.Volume)
hw_forecast_iom <- forecast(hw_model_iom, xreg = test_data_iom$Total_Total.Volume, h = 12)
arima_cme <- auto.arima(price_cme, seasonal = FALSE)
arima_forecast_cme <- forecast(arima_cme, h = 12)
arima_imm <- auto.arima(price_imm, seasonal = FALSE)
arima_forecast_imm <- forecast(arima_imm, h = 12)
arima_iom <- auto.arima(price_iom, seasonal = FALSE)
arima_forecast_iom <- forecast(arima_iom, h = 12)
# price_cme, price_imm, price_iom already has frequency of 12
sarima_cme <- auto.arima(price_cme, seasonal = TRUE)
sarima_forecast_cme <- forecast(sarima_cme, h = 12)
sarima_imm <- auto.arima(price_imm, seasonal = TRUE)
sarima_forecast_imm <- forecast(sarima_imm, h = 12)
sarima_iom <- auto.arima(price_iom, seasonal = TRUE)
sarima_forecast_iom <- forecast(sarima_iom, h = 12)
# Looking at ACF
acf(price_cme)
acf(price_imm)
acf(price_iom)
The autocorrelations show trend, so arfima is applicable.
library(fracdiff)
d_cme <- fracdiff(price_cme)
st_cme <- diffseries(price_cme, d_cme$d)
arfima_cme <- auto.arima(st_cme, xreg = train_data_cme$Total_Total.Volume)
arfima_forecast_cme <- forecast(arfima_cme, xreg = test_data_cme$Total_Total.Volume, h = 12)
d_imm <- fracdiff(price_imm)
st_imm <- diffseries(price_imm, d_imm$d)
arfima_imm <- auto.arima(st_imm, xreg = train_data_imm$Total_Total.Volume)
arfima_forecast_imm <- forecast(arfima_imm, xreg = test_data_imm$Total_Total.Volume, h = 12)
d_iom <- fracdiff(price_iom)
st_iom <- diffseries(price_iom, d_iom$d)
arfima_iom <- auto.arima(st_iom, xreg = train_data_iom$Total_Total.Volume)
arfima_forecast_iom <- forecast(arfima_iom, xreg = test_data_iom$Total_Total.Volume, h = 12)
#install.packages("fGarch")
library(fGarch)
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")
num_train_data_cme <- scale(select_if(train_data_cme, is.numeric))
num_train_data_imm <- scale(select_if(train_data_imm, is.numeric))
num_train_data_iom <- scale(select_if(train_data_iom, is.numeric))
num_test_data_cme <- scale(select_if(test_data_cme, is.numeric))
num_test_data_imm <- scale(select_if(test_data_imm, is.numeric))
num_test_data_iom <- scale(select_if(test_data_iom, is.numeric))
garch_model_cme <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_cme, trace = FALSE)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
garch_forecast_cme <- predict(garch_model_cme, data = num_test_data_cme, n.ahead = 12)
garch_model_imm <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_imm, trace = FALSE)
garch_forecast_imm <- predict(garch_model_imm, data = num_test_data_imm, n.ahead = 12)
garch_model_iom <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_iom, trace = FALSE)
garch_forecast_iom <- predict(garch_model_iom, data = num_test_data_iom, n.ahead = 12)
cme_rmse <- data.frame(Commodity = "CME",
LM = accuracy(lm_forecast_cme, test_data_cme$Price)[[2]],
ARMA = accuracy(arma_forecast_cme, test_data_cme$Price)[[2]],
HW = accuracy(hw_forecast_cme, test_data_cme$Price)[[2]],
ARIMA = accuracy(arima_forecast_cme, test_data_cme$Price)[[2]],
SARIMA = accuracy(sarima_forecast_cme, test_data_cme$Price)[[2]],
ARFIMA = accuracy(arfima_forecast_cme, test_data_cme$Price)[[2]],
GARCH = accuracy(garch_forecast_cme$meanForecast, num_test_data_cme[,1])[[2]])
imm_rmse <- data.frame(Commodity = "IMM",
LM = accuracy(lm_forecast_imm, test_data_imm$Price)[[2]],
ARMA = accuracy(arma_forecast_imm, test_data_imm$Price)[[2]],
HW = accuracy(hw_forecast_imm, test_data_imm$Price)[[2]],
ARIMA = accuracy(arima_forecast_imm, test_data_imm$Price)[[2]],
SARIMA = accuracy(sarima_forecast_imm, test_data_imm$Price)[[2]],
ARFIMA = accuracy(arfima_forecast_imm, test_data_imm$Price)[[2]],
GARCH = accuracy(garch_forecast_imm$meanForecast, num_test_data_imm[,1])[[2]])
iom_rmse <- data.frame(Commodity = "IOM",
LM = accuracy(lm_forecast_iom, test_data_iom$Price)[[2]],
ARMA = accuracy(arma_forecast_iom, test_data_iom$Price)[[2]],
HW = accuracy(hw_forecast_iom, test_data_iom$Price)[[2]],
ARIMA = accuracy(arima_forecast_iom, test_data_iom$Price)[[2]],
SARIMA = accuracy(sarima_forecast_iom, test_data_iom$Price)[[2]],
ARFIMA = accuracy(arfima_forecast_iom, test_data_iom$Price)[[2]],
GARCH = accuracy(garch_forecast_iom$meanForecast, num_test_data_iom[,1])[[2]])
combined_rmse <- rbind(cme_rmse, imm_rmse, iom_rmse)
combined_rmse$Best_Model <- apply(combined_rmse[, -1], 1, function(x) {
colnames(combined_rmse[, -1])[which.min(abs(x))]
})
print(combined_rmse)
## Commodity LM ARMA HW ARIMA SARIMA ARFIMA
## 1 CME 271928.7 -241806.2206 -192.5793 31291.667 31291.667 512777.13
## 2 IMM 242626.0 25469.1639 33720.7608 34281.914 34281.914 268904.23
## 3 IOM 116316.6 -900.4912 14852.5595 -1071.768 -4836.629 89559.59
## GARCH Best_Model
## 1 1.064806 GARCH
## 2 1.242325 GARCH
## 3 1.606774 GARCH
Best Model: - GARCH for all
sts <- read.csv(paste(dataPath,'sts_data_2023.csv',sep = '/'), header=TRUE)
sts <- na.omit(sts)
sts_1990 <- sts %>%
filter(as.numeric(format(as.Date(Date, "%m/%d/%y"), "%y")) > 23) %>%
mutate(Date = format(as.Date(Date, "%m/%d/%y"), "19%y-%m-%d"))
sts_2000 <- sts %>%
filter(as.numeric(format(as.Date(Date, "%m/%d/%y"), "%y")) <= 23) %>%
mutate(Date = format(as.Date(Date, "%m/%d/%y"), "20%y-%m-%d"))
sts <- rbind(sts_1990, sts_2000)
sts <- xts(sts[,-1], order.by = as.Date(sts$Date))
#install.packages("bsts")
library(bsts)
## Loading required package: BoomSpikeSlab
## Loading required package: Boom
##
## Attaching package: 'Boom'
## The following object is masked from 'package:stats':
##
## rWishart
##
## Attaching package: 'BoomSpikeSlab'
## The following object is masked from 'package:stats':
##
## knots
##
## Attaching package: 'bsts'
## The following object is masked from 'package:BoomSpikeSlab':
##
## SuggestBurn
data_fit1 <- window(sts, end = "1979-06-30")
data_fit2 <- window(sts,end = "2004-06-30")
data_fit3 <- window(sts, end = "2022-12-31")
new_data1 <- head(window(sts, start = "1979-09-30"), 5)
new_data1_CL <- new_data1[, -3]
new_data1_GC <- new_data1[, -1]
new_data2 <- head(window(sts, start = "2004-09-30"), 5)
new_data2_CL <- new_data1[, -3]
new_data2_GC <- new_data1[, -1]
new_data3 <- head(window(sts, start = "2023-03-31"), 1)
new_data_3CL <- new_data1[, -3]
new_data3_GC <- new_data1[, -1]
# FIT 1
ss <- AddLocalLevel(list(), data_fit1$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit1$Mean_CL ~ data_fit1$Mean_GC + data_fit1$Mean_CPI + data_fit1$Mean_SPX)
model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit1, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:43 2024 =-=-=-=-=
plot(model_dynamic)
pred_CL_dynamic_1 <- predict.bsts(model_dynamic, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_dynamic_1)
# FIT 2
ss <- AddLocalLevel(list(), data_fit2$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit2$Mean_CL ~ data_fit2$Mean_GC + data_fit2$Mean_CPI + data_fit2$Mean_SPX)
model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit2, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:44 2024 =-=-=-=-=
plot(model_dynamic)
pred_CL_dynamic_2 <- predict.bsts(model_dynamic, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_dynamic_2)
# FIT 3
ss <- AddLocalLevel(list(), data_fit3$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit3$Mean_CL ~ data_fit3$Mean_GC + data_fit3$Mean_CPI + data_fit3$Mean_SPX)
model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit3, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_dynamic)
pred_CL_dynamic_3 <- predict.bsts(model_dynamic, newdata = new_data3, horizon = 1, quantiles = c(.025, .975))
plot(pred_CL_dynamic_3)
# FIT 1
ss <- AddLocalLinearTrend(list(), data_fit1$Mean_CL)
model_trend <- bsts(data_fit1$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_trend)
pred_CL_trend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_1)
# FIT 2
ss <- AddLocalLinearTrend(list(), data_fit2$Mean_CL)
model_trend <- bsts(data_fit2$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_trend)
pred_CL_trend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_2)
# FIT 3
ss <- AddLocalLinearTrend(list(), data_fit3$Mean_CL)
model_trend <- bsts(data_fit3$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=
plot(model_trend)
pred_CL_trend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_3)
# FIT 1
ss <- AddSemilocalLinearTrend(list(), data_fit1$Mean_CL)
model_trend <- bsts(data_fit1$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=
plot(model_trend)
pred_CL_semitrend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_1)
# FIT 2
ss <- AddSemilocalLinearTrend(list(), data_fit2$Mean_CL)
model_trend <- bsts(data_fit2$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=
plot(model_trend)
pred_CL_semitrend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_2)
# FIT 3
ss <- AddSemilocalLinearTrend(list(), data_fit3$Mean_CL)
model_trend <- bsts(data_fit3$Mean_CL, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=
plot(model_trend)
pred_CL_semitrend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_3)
# FIT 1
ss <- AddLocalLevel(list(), data_fit1$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit1$Mean_GC ~ data_fit1$Mean_CL + data_fit1$Mean_CPI + data_fit1$Mean_SPX)
model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit1, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:47 2024 =-=-=-=-=
plot(model_dynamic)
pred_GC_dynamic_1 <- predict.bsts(model_dynamic, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_dynamic_1)
# FIT 2
ss <- AddLocalLevel(list(), data_fit2$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit2$Mean_GC ~ data_fit2$Mean_CL + data_fit2$Mean_CPI + data_fit2$Mean_SPX)
model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit2, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:47 2024 =-=-=-=-=
plot(model_dynamic)
pred_GC_dynamic_2 <- predict.bsts(model_dynamic, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_dynamic_2)
# FIT 3
ss <- AddLocalLevel(list(), data_fit3$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit3$Mean_GC ~ data_fit3$Mean_CL + data_fit3$Mean_CPI + data_fit3$Mean_SPX)
model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit3, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=
plot(model_dynamic)
pred_GC_dynamic_3 <- predict.bsts(model_dynamic, newdata = new_data3, horizon = 1, quantiles = c(.025, .975))
plot(pred_GC_dynamic_3)
# FIT 1
ss <- AddLocalLinearTrend(list(), data_fit1$Mean_GC)
model_trend <- bsts(data_fit1$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=
plot(model_trend)
pred_GC_trend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_1)
# FIT 2
ss <- AddLocalLinearTrend(list(), data_fit2$Mean_GC)
model_trend <- bsts(data_fit2$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=
plot(model_trend)
pred_GC_trend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_2)
# FIT 3
ss <- AddLocalLinearTrend(list(), data_fit3$Mean_GC)
model_trend <- bsts(data_fit3$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)
pred_GC_trend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_3)
# FIT 1
ss <- AddSemilocalLinearTrend(list(), data_fit1$Mean_GC)
model_trend <- bsts(data_fit1$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)
pred_GC_semitrend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_1)
# FIT 2
ss <- AddSemilocalLinearTrend(list(), data_fit2$Mean_GC)
model_trend <- bsts(data_fit2$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)
pred_GC_semitrend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_2)
# FIT 3
ss <- AddSemilocalLinearTrend(list(), data_fit3$Mean_GC)
model_trend <- bsts(data_fit3$Mean_GC, state.specification = ss, niter = 200)
## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=
plot(model_trend)
pred_GC_semitrend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_3)
CL_fit1_rmse <- data.frame(Commodity = "CL Fit 1",
DynamicRegression = accuracy(pred_CL_dynamic_1$mean, new_data1$Mean_CL)[[2]],
LinearTrend = accuracy(pred_CL_trend_1$mean, new_data1$Mean_CL)[[2]],
SemiLocalTrend = accuracy(pred_CL_semitrend_1$mean, new_data1$Mean_CL)[[2]])
CL_fit2_rmse <- data.frame(Commodity = "CL Fit 2",
DynamicRegression = accuracy(pred_CL_dynamic_2$mean, new_data2$Mean_CL)[[2]],
LinearTrend = accuracy(pred_CL_trend_2$mean, new_data2$Mean_CL)[[2]],
SemiLocalTrend = accuracy(pred_CL_semitrend_1$mean, new_data2$Mean_CL)[[2]])
CL_fit3_rmse <- data.frame(Commodity = "CL Fit 3",
DynamicRegression = accuracy(pred_CL_dynamic_3$mean, new_data3$Mean_CL)[[2]],
LinearTrend = accuracy(pred_CL_trend_3$mean, new_data3$Mean_CL)[[2]],
SemiLocalTrend = accuracy(pred_CL_semitrend_3$mean, new_data3$Mean_CL)[[2]])
GC_fit1_rmse <- data.frame(Commodity = "GC Fit 1",
DynamicRegression = accuracy(pred_GC_dynamic_1$mean, new_data1$Mean_GC)[[2]],
LinearTrend = accuracy(pred_GC_trend_1$mean, new_data1$Mean_GC)[[2]],
SemiLocalTrend = accuracy(pred_GC_semitrend_1$mean, new_data1$Mean_GC)[[2]])
GC_fit2_rmse <- data.frame(Commodity = "GC Fit 2",
DynamicRegression = accuracy(pred_GC_dynamic_2$mean, new_data2$Mean_GC)[[2]],
LinearTrend = accuracy(pred_GC_trend_2$mean, new_data2$Mean_GC)[[2]],
SemiLocalTrend = accuracy(pred_GC_semitrend_1$mean, new_data2$Mean_GC)[[2]])
GC_fit3_rmse <- data.frame(Commodity = "GC Fit 3",
DynamicRegression = accuracy(pred_GC_dynamic_3$mean, new_data3$Mean_GC)[[2]],
LinearTrend = accuracy(pred_GC_trend_3$mean, new_data3$Mean_GC)[[2]],
SemiLocalTrend = accuracy(pred_GC_semitrend_3$mean, new_data3$Mean_GC)[[2]])
combined_rmse <- rbind(CL_fit1_rmse, CL_fit2_rmse, CL_fit3_rmse, GC_fit1_rmse, GC_fit2_rmse, GC_fit3_rmse)
combined_rmse$Best_Model <- apply(combined_rmse[, -1], 1, function(x) {
colnames(combined_rmse[, -1])[which.min(abs(x))]
})
print(combined_rmse)
## Commodity DynamicRegression LinearTrend SemiLocalTrend Best_Model
## 1 CL Fit 1 30.39734 6.255817 12.663403 LinearTrend
## 2 CL Fit 2 38.35451 3.689569 28.469102 LinearTrend
## 3 CL Fit 3 78.29768 2.932757 5.081877 LinearTrend
## 4 GC Fit 1 458.76358 249.223872 266.072375 LinearTrend
## 5 GC Fit 2 347.42962 11.774606 134.081005 LinearTrend
## 6 GC Fit 3 1531.99188 195.287713 204.353104 LinearTrend
QUESTIONS Local linear trend works best with everything except fit 3 of GC, which is best modeled using Semi Local Linear Trend state-specification. Local Linear Trend is best overall.
Lags are not useful in predicting CL and GC. All coefficients are 0 meaning no meaningful linear relationship between variables.